#pragma once
#include "feDefs.h"
#include "feSmartPtr.h"
#include <math.h>

class feDenseMatrix;
class feVector;

//-------------------------------------------------------------------------------------------------------------------------------------
class feMatrix
{
public:
	virtual int		GetSize() const = 0;
	virtual feReal*	_IJ(int iI, int iJ) = 0;
	virtual bool	LuDecomposition() = 0;
	virtual bool	ForwardTrace(feVector *iF, feVector *oX) = 0;
	virtual bool	BackwardTrace(feVector *iF, feVector *oX) = 0;
	virtual void	AddM(feMatrix *iMtoAdd, double iAlpha=1) {  }
	virtual void	MultMv(feVector *iV, feVector *oV) = 0;
	virtual void	MultMAlpha(feReal iAlpha) = 0;
	virtual void	AddLocalToGlobal(feDenseMatrix *iM, int *iMRemap, double iMultCoeff) = 0;
	virtual void	Zero() = 0;
	virtual ~feMatrix() {};
};


//-------------------------------------------------------------------------------------------------------------------------------------
class feDenseMatrix : public feMatrix
{
	feReal **A;
	int		 N;


private:
	void Free();


public:
	feDenseMatrix(int iN);
	feDenseMatrix() { A = 0; N = 0;}

	virtual int		GetSize() const		{ return N; }
	virtual double*	_IJ(int iI, int iJ)	{ return &A[iI][iJ]; }
	virtual bool	LuDecomposition();
	virtual bool	ForwardTrace(feVector *iF, feVector *oX);
	virtual bool	BackwardTrace(feVector *iF, feVector *oX);
	virtual void	MultMv(feVector *iV, feVector *oV);
	virtual void	MultMAlpha(feReal iAlpha);
	virtual void	AddLocalToGlobal(feDenseMatrix *iM, int *iMRemap, double iMultCoeff) {}
	virtual void	Zero();
			void	Fill(feReal **iFrom);
			void	Transpose(feDenseMatrix *oMatrix);
			void	Square(feDenseMatrix *oM);
			void	Allocate(int iSize);
	~feDenseMatrix();

};


//-------------------------------------------------------------------------------------------------------------------------------------
class feFiveDiagSymMatrix : public feMatrix
{
	feReal *A[3];
	int		 N;

public:
	feFiveDiagSymMatrix();
	feFiveDiagSymMatrix(int iN);
	feFiveDiagSymMatrix(feFiveDiagSymMatrix& iM);

			void	Allocate(int iN);
	virtual int		GetSize() const			{ return N; }
	virtual feReal*	_IJ(int iI, int iJ);
	virtual void	AddM(feMatrix *iMtoAdd, double iAlpha=1);
	virtual bool	LuDecomposition();		
	virtual bool	ForwardTrace(feVector *iF, feVector *oX);
	virtual bool	BackwardTrace(feVector *iF, feVector *oX);
	virtual void	MultMv(feVector *iV, feVector *oV);
	virtual void	MultMAlpha(feReal iAlpha);
	virtual void	AddLocalToGlobal(feDenseMatrix *iM, int *iMRemap, double iMultCoeff);
	virtual void	Zero();
			void	Free();

	~feFiveDiagSymMatrix();
};


//-------------------------------------------------------------------------------------------------------------------------------------
class feThreeDiagSymMatrix : public feMatrix
{
	feReal *A[2];
	int		 N;


public:
	feThreeDiagSymMatrix();
	feThreeDiagSymMatrix(int iN);
	feThreeDiagSymMatrix(feThreeDiagSymMatrix& iM);

	void			Allocate(int iN);
	virtual int		GetSize() const			{ return N; }
	virtual feReal*	_IJ(int iI, int iJ);
	virtual void	AddM(feMatrix *iMtoAdd, double iAlpha=1);
	virtual bool	LuDecomposition();		
	virtual bool	ForwardTrace(feVector *iF, feVector *oX);
	virtual bool	BackwardTrace(feVector *iF, feVector *oX);
	virtual void	MultMv(feVector *iV, feVector *oV);
	virtual void	MultMAlpha(feReal iAlpha);
	virtual void	AddLocalToGlobal(feDenseMatrix *iM, int *iMRemap, double iMultCoeff);
	virtual void	Zero();
	void			Free();

	~feThreeDiagSymMatrix();
};


